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Effects of non-Gaussian a— stable Levy noise on the Gompertz tumor growth model are quan- 
tified by considering the mean exit time and escape probability of the cancer cell density from 
inside a safe or benign domain. The mean exit time and escape probability problems are formu- 
lated in a differential-integral equation with a fractional Laplacian operator. Numerical simula- 
tions are conducted to evaluate how the mean exit time and escape probability vary or bifurcates 
when a changes. Some bifurcation phenomena are observed and their impacts are discussed. 
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1. Introduction 

Mathematical models have been proposed to describe the evolution of numerous biological phenomena such 
as cancer growth. The deterministic Gompertz growth model for a tumor density growth is given by the 
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following differential equation 



— — = (A — Bmx]x, 
dt ^ ' ' 



(1) 



where x(t) describes the density of cancer cells at time t, parameters A and B denote growth and decay 
rates respectively. 

This is a deterministic model, nevertheless. As there is a discrepancy between clinical data and basic 
theory due to subtle environmental fluctuations, we need a stochastic part to describe the fluctuation. 
A stochastic model proposed as a fun ctiona l Fokker-Planck equation in consideration of both fission and 
mortality has b een considered in Lol. 12007 ] and a diffusion proc ess with modified infinitesimal mean has 
been studied in [Gutierrez-jaimez et al\.\im% . In [LoI. l2009tl and lAlbano h Giornd . [200911 a multiplicative 
noise of Brownian motion has been considered and in [Albano fc Giornol . l200()f a dete rministic parameter 
has been modified by a group of random variables to depict the stochastic fiuct^ation OB compared 
the effec t of different therapies by giving their probability density functions(pdfs) and Albano &: Giornol . 
I2OO9], [ Albano fc Giornol . 120061 ] both deal with the first exit time problem. 

In the present paper we consider the Gompertz growth model driven by the non-Gaussian a— stable 
noise with jumps and quantify the dynamics in terms of mean exit time and escape probability. Specifically, 
we consider the evolution of the tumor growth in the density range that cannot be diagnosed, and compute 
the likelihood and mean time for a tumor to become diagnosable or malignant. We examine how this 
depends on a values. To this end, we investigate the differential equation driven by symmetric a— stable 



Levy process with characteristics (0, 1, z^) with jump measure v(dy) = dy/\y\^ 



dX = [{A -BlnX)X]dt + dL 



a 

t 5 



A(0) 



X. 



(2) 



where x is the initial tumor density. 

The paper is arranged as follows. We present a brief introduction about Levy process in section 2 and 
then a stochastic tumor growth model in section 3. Numerical experiments are conducted in section 4. 



2. Levy process 

A Levy process (Lt, t > 0) is a stochastic process defined on a probability space (0, T, P) with the following 
properties: 

(i) Lo = (a.s.); 

(ii) Lt has independent and stationary increments, i.e. for each n and < ti < t2 < ■ ■ ■ < tn+i < 00 the 
random variables {Lt^^^ — Lt- ,1 < j < n) are independent and each Ltj_^_^ — Lt - and Lt-j^^-tj are equal 
in law both satisfy a non-Gaussian distribution. 

(iii) Lt has stochastically continuous, i.e. for all c > 0, s > 

\i-mP{\Lt - Ls\ > c) = 0. 

(iv) The paths of Lt are P— almost surely right continuous with left limits. 

For the characteristic exponent ^t(A) = logE(e*'^^*), with the fact that Lt has independent and 
stationary increment we get E(e*''^^') = e**^^"**^ =: e**^'^\ therefore in the next paragraph we refer to ^'(A) 
as a characteristic exponent of the Levy process . 

Levy-Khintchine formula |Applebaum . 20091 ; IKyprianou . 20061 ] for Levy process: Suppose 

that a € M, o" > and v is a. measure concentrated on M \ {0} such that A x^) fj,{dx) < 00. From this 
triple (a, a, v) defined for each A € M, 

^(A) = ia\ - ^(j^A^ - / (1 - e^^'^ + «Ax/|^|<i) v{dx). 

Then there exits a probability space (17,7-", P) on which a Levy process is defined having characteristic 
exponent ^ . 
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It can be seen that after some reorganization the characteristic exponent of Levy process can be 
rewritten as 

-\ [ {l-e'^'-" + iXx)u{dx)} 

*- Jo<\x\<l J 



'0<|a;|<l 

^(1) _ ^(2) _ ^(3) 



where ^'^^^ are respectively the characteristic exponent of xj:^^ = —aBt + at, t > and xj:"^^ 

X^i^i^i) t > 0, here {Nt : t > 0} is a Poisson process with rate z^(M \ (—1, 1)) and {Ci ■ t > 0} are i.i.d 

(2) 

with distribution ij{dx)/u(M. \ (—1, 1)) concentrated on {x : \x\ > 1} (if \ (—1, 1)) = process X^ is 
identicaUy zero). 

We have a result from [Kvprianou . 2006 ] to say the composition of a Levy process. 



Levy-Ito decomposition: Given any o € M, cr > and measure v concentrated on M\ {0} satisfying 
/r\{o}(I^I^ a l)i'{dx) < oo, there exists a probability space on which three independent Levy processes 
exit, where X^^) is a linear Brownian motion with drift with characteristic exponent ' , 
X^'^'^ is a compound Poisson process with characteristic exponent ^^"^^ and X^^^ is a square integrable 
martingale with an almost surely countable number of jumps on each finite time interval which are of 
magnitude less than unity and with characteristic exponent By taking X = X^^^ — X^'^^ — X^^^ we 

have that there exist a probability space on which a Levy process is defined with characteristic exponent 

*(A) = iaX - Ja^A^ - / (1 - e^^^ + iAx/(|^|<i)) fi{dx) 
2 Jr 



for A € 



A random variable X is said to be stable if for all n > 0, the distribution of X and S'„ 



\-Xn —an 



0"n 

are equivalent, whe re X^, - ■ ■ , Xn are a sequence o f i.i.d random variable, a„ € R and (t„ > for n € N. It 
is well known from Applebaum . 20091 : Feller . 1971 1 that the only possible choice is o"„ = n^/" for a E (0, 2] 



so we refer to a as the index. 

For a E (0, 1) U (1,2) stable variables have characteristic exponents of the form 

^{X) = iaA-c"|A|"(l-i/3tan— s^nA), 

where (3 G [—1, 1], a G M, c > and the sign function sgnX = I(x>o) ~ ^{x<o)'^ 
for a = 1, the characteristic exponents have the form 

2 

qi(X) =mA-c|A|(l + z/3-sc/nAlog|A|), 

TT 



where f3 G [-1, 1], a G M and c > 0; 
for the special case a = 2, 



and the measure 

i^{dx) 



*(A) = iaX - ^CJ^A^; 



-dx for X G (0, oo); 



jzi^frr^dx for xG(-oo,0), 



where c = ci + C2, ci, C2 > and /3 = if a G (0, 1) U (1, 2). The case /3 = 0, i.e. ci = C2, is called 

symmetric a— stable Levy process. 
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The (infinitesimal) generator A of a Levy process X{t) is defined as Aip = lim '^"'^ where {Ttip){x) 



t io 



E,x{f{X{t)) for any ip € -D^(the domain of generator A). It has the expression [Applebauml . I2nn9l | 

{Aip){x) = a(p'{x) + ]-a^ip"{x) + / [ip{x + y) - ip{x) - yv?'(x)/|y[<i] v{dy). 



3. Tumor growth under a— stable Levy noise 

In this section we consider the stochastic model driven by symmetric a— stable Levy process, i.e. a Levy 
process Lj, in which each Lt is a stable random variable and here for simplicity we take ci = C2 = 1. We refer 
to the monoclonal benign parathyroid tumor and take A = 6.46 year~^, B = 0.314 year~^, corresponding 
to a parathyroid tumor mean age of 19.6 years. Here we investigate in the range of the tumor density that 
cannot be diagnosed, we take the two boundary points Si = 1 and 52 = 1.074 x 10^ (the tumor cell density 
X = 1.074 X 10^ is the smallest diagnoseable mass), how the tumor density progresses. When a tumor's 
density is bigger than 1.074 x 10^, we say it becomes malignant, while if its density is smaller than 1, it is 
not diagnosable. 

Based on the deterministic Gompertz model, the corresponding stochastic model is given by 



dX = [{A- B\nX)X]dt + dL^, X(0) = x, 

where L" is with the characteristics (0,1, z^). Here the jump measure i'{dy) = dy/\y\^^"'. 
(Applebaum . 20091 ] that the infinitesimal generator A for Lf is 



{Af){x) 



f{x + y)-f{x)-yf'{x)Iu^i 



\l+a 



/M\{o} \y\ 

Due to the Levy-Ito decomposition the generator for the process X{t) in 



dy. 



is then 



Af = {A-B In x)x fix) + ^f"{x) + 



fix + y)- fix) - yf'ix)I\y\^i 



\l+a 



dy. 



(3) 

It is known 

(4) 

(5) 



In this symmetric a— stable Levy motion case, th e integral operato r is related to the frac- 
tional Laplacian operator lAlbeverrio. Rudiger fc Wul . l200d : iGuan fc Mai . l2005l : IChen. Kim Soii3 . l20ld : 
Caffarelh fc Silvestrd . l2007l ] . Namelv. for a F (0.21 

/ [uix + y) - n(x)] i^aidy) = Ca (-A)tu(x), (6) 
Jr\{o} 

where 

(cosy - 1) i^aidy) < 0. 



Note that when a = 2, L2 = A. 

We consider two issues: Mean exit time and escape probability. We quantify the stochastic dynamics 
using these two tools. Especially, we examine possible bifurcation phenomena when the stability parameter 
a various in (0, 2). 

The mean exit time 

uix) = E mf{t > 0, Xtiuj, x) ^ D}, 

from a bounded interval D, is the first time that the tumor density Xit), initially at value x in D, gets 
out of D. The mean exit time ujx), from the diagnos eable range D = (1, 1.074 x 10^), satisfies an integral- 
differential equation Brannan. Duan &: Ervin . 19991 ] 



Auix) 
uix) 



-1, x£ D 
0, x^D. 



(1,1.074 X 10^), 



(7) 
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i.e. 

{A-B\nx)xu'{x) + -u"{x)+ / — , ' ^ ' '^'^ dy = -1, xeD, 

u{x) = 0, x^D. (8) 

The mean exit time u{x) is the the time when the tumor (initiahy with density x) becomes not diagnoseable 
(when x{t) exits from the left boundary point x = 1), or becomes mahgnant (when x{t) exits from the 
right boundary point x = b) and thus may need medical intervention. 

The escape probability p{x), through the right boundary point, is the likelihood that the tumor density 
X{t), init ially at value x in -D, first escape s from D through x = 1.074 x 10^. The p{x) satisfies the following 
equation jBrannan. Duan &: Erviiil . Il999l | 

Ap{x) = 0, X e D = (1,1.074 X 10^), 
p{l) = 0, p(1.074 X 10*^) = 1. (9) 

In the next section we simulate both the mean exit time u{x) and escape probability p{x), when the 
stability parameter a varies in (0,2). 



4. Numerical experiments 

We now discretize the equation ([8]) for the mean exit time u{x). The discretization of ([9]) for the escape 
probability p{x) is similar. 

For the numerical scheme as the right boundary point of D is too huge, we need a transformation to 
change it to a smaller value. Owing to the fact that u{x) = for x ^ D we can take a transformation 
X = Inx, relevantly y = ln(x + y) — Inx, nonetheless the integration for y is on M \ {0} {u{x + y) = for 
X + y < < 1). Omitting the tilde we get 

{A - Bx)u'{x) + ^e-2^(u"(x) - n'(x)) 

1 r uix + y) — uix) — ie^ — l)u'{x)Is\f.yT\^p-x\ 
+ 4: , , ^' K ydy=-1, xGDi = (0,18.4921), 



u(x) = 0, x^Di. (10) 

where the third term in the integral is 

\^y_l\i+a Jmsnl |e?'-ir+" ■Ir\{0} W] 



as /]j\|o} ^l^jiVlf* dy always vanishes for any 5 > 0, so does the third term in the integral for any fixed x, 
therefore the integral-differential equation can be changed into 

^e-2^'u"(x) + [{A - Bx) - ie-2^]n'(x) 

+ + " " -l>'^-)h\ey-l\<S) ^^^^^ 

u(x) = 0, xiDi. (11) 

We then take another transformation x = | — 1 meanwhile y = | where h = 18.4921/2 to change 
Di = (0, 18.4921) into (-1, 1), we obtain 

_Le-2K-+i)n"(x) +[A- Bb{x + 1) - le'^^^^+i)] ^n'(x) 

I 1 [ fe^(x + y)- bujx) - {e'y - l)u'{x)I{^,.y^,^^s} J,, _ 

n(x) = 0, X ^ D2. (12) 



,ab(x+l) 1^^^^^ lefts, -l\l+c^ 
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The numerical scheme here fohows jLi. Gao &: Duaiil . 1201 A different scheme was also pro- 
posed in (Chen. Duan. Li fc Zhand . 12011 ]. We write Jj^ = Jl^^ + + ^ ~ 
_ 1^ 1 _ e-f'(i+^)} to get 



mm 



1 g-26(x+i)^//(^) + - Bb{x + 1) - ie-2''(-+i) 
2o^ L 2 



-u (x) 



+ 



+ 



u{x) 



1 



1 



1 



1 



i-^' + y) - 6n(x) - {e^ - l)u' {x)I{\^by_i\^s} j,y 



l-x 



furthermore, for x > 
1 



262' 



+ 



-2b{x+l)„fl 



u{x) 



u"{x) + 



1 



1 



(13) 



A - Bh{x + 1) - ^e-2''(^+i) 



-u'{x) 



1 



1 



+ 



Q,go6(a;+l) [ ||gb(l-x) _ _ g-b(l+a;)^a 

TTT / 

l+x 



+ 



1 



,ab{x+l) 



-l+x 



-l-x 



u{x + y) — u{x) 

\gby _ \\l+oi 



be^y dy 



bu{x + y)- bu{x) - {e^ - l)u'{x) ^ ( ^ ^^ 

— ui 1 1N / rrn e ^ ay = —1, x£lJ3 = (— 1,1) 



(14) 



and for X < 



-Le-2f'(-+i)n"(x) + U - Bb{x + 1) - ie-2''(-+i)l i^,'(x) 
2o^ L 2 J 



+ 



+ 



u{x) 



^^ab{x+l) 
1 

pab{x+l) 



1 



1 



1 



(^eH^-x) _ 1)Q (1 _ g-b(l+a;))c 



+ 



1 



^ab{x+l) 



l-x 



l+x 



u{x + y) — u{x) 

\gby _ \\l+oi 



be^y dy 



^+'^' bu{x + y)- bu{x) - {e^y - l)u'(x) 



l-x 



1+a 



e^dy = -l, xeDs 



■1,1) 



(15) 



For simplicity denote /(x) = A — Bb{x + 1) — |e~2^(^+i) and ^(x) = g-2b(a:+i)^ numerical schemes 
generally we divide the interval [—2, 2] into 4J subintervals and define Xj = jh for —2 J < j < 2 J integer, 
where h = 1/J. We denote the numerical solution of u at xj by Uj. We can discretize the two integral- 
differential equations above using central difference for derivatives and "punched-hole" trapezoidal rule 



+ 



+ 



1 



^pab{xj+l) 

h 

pab{xj+l) 



1 



2h 
1 



(^pb{l-Xj) _ i^a _ p-b{l+Xj)'^a 



h 



,ab{xj + l) 



-J+i 

E 

k=-j-j 



„ be^yf'ia 



j+k - Uj) 



J-j 

E 



-1, (16) 



k=-J+j, k=/=0 

where j = 0, 1, 2, • • • , J — 1. Here the summation symbol Yl " means the two end terms are multiplied by 
1 /2 corresponding to trapezoidformula. 

2b^^^ /i2 + ^A^jj 

J-i 



+ 



1 



+ 



^^ab{xj + l) 

h 



1 



1 



2h 
1 



(^pb{l-Xj) _ i^^a (j^ _ ^-b{l+Xj)-^a 



h 



E 



„ be^'y^iUj+k-Uj 



+ pab{xj+l) |g6j/fe - 1|1+" 

k=J+j 



gab{xj+l) 



-J-j 

E 

k=J+j, k^O 



„ be'y^{Uj+k - Uj) - e'^y^je'^y^ - l){Uj+i - U,-i)/2h 

\gbyk _ 



-1, 



(17) 
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where j = —(J — 1), — ( J — 2), • • • , —2, —1. The boundary conditions require that the values of Uj vanish 
if the index \j\ > J. But here for the specific problem for first-order derivative we use forward difference 
or backward difference. 

Similarly, we can discretize the equation ([9]) for the escape probability p{x). 

For the viewing convenience, the following figures for the mean exit time u{x) and escape probability 
p{x) are both plotted for x £ (0, 18.4921), instead of the original huge range (1, 1.074 x 10^). 

Figures [iHH] show the mean exit time u{x) for various a values in (0,2). 

We observe that when a is small, the mean exit time u{x) is small and is roughly constant, i.e., tumors 
at all densities are equally unlikely to become non-diagnosable or malignant. As a becomes bigger the 
whole height will grow with the right part grows up slowly while the left part increases rapidly meanwhile 
the width of the left part narrows down. In fact, at about a = 0.2 the height of the two part begins to have 
a separation. When a ~ 0.44 the separation becomes quite obvious and the left part (for density between 
and 10) reaches its highest value, which indicates a bifurcation. This means that for a tumor with density 
between and 10, it takes longest time to become either non-diagnosable or malignant. 

For a between 0.44 and 0.48, u{x) keeps similar shape as for a = 0.44 but the peak diminishes. For 
a between 0.48 and 0.6, u{x) has the similar shape although the peak value decreases. Moreover, for a 
between 0.6 and 2, u{x) keeps almost the same shape; namely, for the Gaussian Brownian motion case, the 
mean exit time is similar to the large values of a € (0.6,2). The dramatic impact of non-Gaussian Levy 
noise occurs for a G (0,0.46). 

Figure [7] shows a three dimensional plot of the above discussed situations. 

In the above simulations, the mean exit time u{x) tells us the time for a tumor to become either non- 
diagnosable (exit from the left boundary point) or become malignant (exit from the right boundary point). 
To distinguish these two situations, especially examine the situation for a tumor to become malignant, let 
us now consider the escape probability, p{x), through the right end point. It indicates the likelihood the 
the tumor is progressing from benign to malignant. We observe that a bifurcation occurs for a ~ 0.4 when 
a tumor at any benign density will highly likely becomes malignant, while for a < 0.4, only tumors with 
high density (near x = 20) are likely to become malignant; see Figures [8] and [9j 

Acknowledgements. We thank Xiaofan Li and Ting Gao for help with numerical schemes, and Xu Sun 
for helpful discussions on the model. 
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Mean Exit Time winen a=0.1 
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Fig. 1. Plot of u{x) given by (|T0|) with a = 0.01 and a = 0.1. 
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Mean Exit Time when a=0.46 
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Fig. 3. Plot of u{x) given by flOj with a = 0.44 and a = 0.46. 
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Mean Exit Time when a=0.48 
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Mean Exit Time when a=0.6 
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Fig. 4. Plot of u{x) given by (fTO)) with a = 0.48 and a = 0.6. 
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Mean Exit Time when a=0.7 
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Fig. 5. Plot ol u{x,t) given by (|10p with a = 0.7 and a — 1. 
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Mean Exit Time when a=1 .8 
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Fig. 6. Plot of u{x,t) given by (|10[) witli a — 1.8 and the Brownian motion 
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Mean Exit Time 




Fig. 7. Mean exit time u v.s. a and x. 
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Fig. 8. 



Escape probability p{x) through the right boundary of (0, 18.4921): 



a — 0.1 and a 



= 0.2. 
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Escape Probablity through the right boundary when a=0.4 
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Escape Probablity through the right boundary when a=0.5 
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Fig. 9. Escape probability p{x) through the right boundary of (0, 18.4921): a — 0.4 and a = 0.5. 



